function ConvStudyMesh1(varargin)
clc
%close all
addpath('../');
addpath('../../');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% program interface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 0
    help();
    return
end
problemType = 0;
echoLevel = 0;
if nargin > 0
    problemType = varargin(1)
end
if nargin > 1
    echoLevel = varargin(2)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create mesh/read mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l = 1;
h = 1;
d = 0.1;
if strcmp(problemType,'simple mesh')
    mesh = fem2.Mesh();
    n1 = mesh.addNode(fem2.Node(0,0));
    n2 = mesh.addNode(fem2.Node((l-d)/2,0));
    n3 = mesh.addNode(fem2.Node((l+d)/2,0));
    n4 = mesh.addNode(fem2.Node(l,0));
    n5 = mesh.addNode(fem2.Node(l,h));
    n6 = mesh.addNode(fem2.Node((l+d)/2,h));
    n7 = mesh.addNode(fem2.Node((l-d)/2,h));
    n8 = mesh.addNode(fem2.Node(0,h));
    mesh.addFace(fem2.Face([n1 n2 n8]));
    mesh.addFace(fem2.Face([n2 n7 n8]));
    mesh.addFace(fem2.Face([n2 n3 n7]));
    mesh.addFace(fem2.Face([n3 n6 n7]));
    mesh.addFace(fem2.Face([n3 n4 n6]));
    mesh.addFace(fem2.Face([n4 n5 n6]));
elseif strcmp(problemType,'mesh2')
    projDir = '/home/hbui/kratos_janosch';
    projName = 'patchtest-conv-study.gid';
    fn = [projDir '/' projName '/' 'mesh2.msh'];
    mg = fem2.gid.GidMeshGenerator(fn);
    mesh = mg.getMesh2d();
elseif strcmp(problemType,'mesh1')
    projDir = '/home/hbui/kratos_janosch';
    projName = 'patchtest-conv-study.gid';
    fn = [projDir '/' projName '/' 'mesh1.msh'];
    mg = fem2.gid.GidMeshGenerator(fn);
    mesh = mg.getMesh2d();
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% config for mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mesh.setNumberOfFields(2);
%set id for nodes
for i=1:mesh.countNodes()
    mesh.getNode(i-1).setId(i);
end
%set id for faces
for i=1:mesh.countFaces()
    mesh.getFace(i-1).setId(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
du = 1;
lambda = linspace(0,1,100);
meshUtil = fem2.MeshUtilities(mesh);
c1 = fem2.Constraint([0 0]); %fix x y
c2 = fem2.Constraint([0 1]); %fix x
c3 = fem2.Constraint([0 0]); %fix y, prescribe x
c3.setValue(0,du);
c4 = fem2.Constraint([1 0]); %fix y
c5 = fem2.Constraint([0 1]); %prescribe x
c5.setValue(0,du);
leftedge = meshUtil.seekNodesOnSurface(1,0,0);
rightedge = meshUtil.seekNodesOnSurface(1,0,-l);
botedge = meshUtil.seekNodesOnSurface(0,1,0);
mesh.setConstraint(c2,leftedge);
mesh.setConstraint(c4,botedge);
meshUtil.seekNode([0 0]).setConstraint(c1);
mesh.setConstraint(c5,rightedge);
meshUtil.seekNode([l 0]).setConstraint(c3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% view the mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(echoLevel,'see mesh')
    plotMesh(mesh);
    return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element declaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = fem2.Model(mesh);
nf = mesh.countFaces();
E = 2;
nu = 0.3;
H = -1;
ft = 0.91;
C = C_plane_stress(E,nu);
for i=1:nf
    e = fem2.MatlabStructuralElement(mesh.getFace(i-1));
    M.addElement(e);
    %create struct to contain element state variables
    es = createElementContainer(C,H,ft);
    elements(i) = es;
end
%deactivate and activate elements
for i=1:nf
    c = computeFaceCenter(mesh.getFace(i-1));
    if (c(1)>=(l-d)/2 && c(1)<=(l+d)/2 ...
            && c(2)>=0 && c(2)<=h)
        elements(i).isActivated = 1;
    else
        elements(i).isActivated = 0;
    end
end

for i=1:leftedge.size()
    leftedgenodes(i) = leftedge.get(i-1).getId();
end

postProcessInfo.infoList = {'resistant force','crack opening'};
postProcessInfo.infoParam = {leftedgenodes,[]};

[elements,output] = solving_strategy_static_2(elements,M,...
    lambda,50,1e-6,postProcessInfo);

figure
res_node_1(:,:) = -output.resistant_force(:,1,:)
delta = lambda*du;
plot(delta,res_node_1);

plotMeshWithCrack(M,elements,{});

plotDeformedMesh(mesh,1);

rmpath('../');
rmpath('../../');
end

function help()
disp([...
    'do you see something?' ...
    ]);
end
